Lattice splitting under intermittent flows 
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We study the splitting ol regular square lattices subject to stochastic intermittent flows. Various 
flow patterns are produced by different groupings of the nodes, based on their random alternation 
between two possible states. The resulting flows on the lattices decrease with the number of groups 
according to a power law. By Monte Carlo simulations we reveal how the time span until the occur- 
rence of a splitting depends on the flow patterns. Increasing the flow fluctuation frequency shortens 
this time span which reaches a minimum before rising again due to inertia effects incorporated in 
the model. The size of the largest connected component after the splitting is rather independent 
of the flow fluctuation frequency but slightly decreases with the link capacities. Our findings carry 
important implications for real-world networks, such as electric power grids with a large share of 
renewable intermittent energy sources. 

PACS numbers: 89.75.Hc, 02.50.Ey, 05.10.-a 



I. INTRODUCTION 

Assessing the robustness of networks against failures 
of nodes and links is an essential research topic across 
many scientific disciplines. Examples range from the 
extinction of species in food webs and malfunctions 
in protein networks to the vulnerability of the World 
Wide Web and cascading failures in electric power grids. 
In the last decade, substantial new insights have been 
gained through the application of methods from statis- 
tical physics PHJ]. Random failures as well as targeted 
attacks have been addressed by first studying static prop- 
erties such as different network topologies [5[. Later on, 
load redistribution models have been introduced to bet- 
ter represent networks supporting the flow of a physi- 
cal quantity. For example, the load of a node has been 
defined by its betweenness centrality @, by the total 
number of efficient paths passing through it 0, or en- 
riched with stochastic flux fluctuations Q- While these 
approaches model the failure propagation in a static man- 
ner, the dynamic flow properties have just recently been 
taken into account 

The contribution of this paper is to investigate the 
impact of stochastic intermittent flow patterns on the 
potential occurrence of cascading link failures, eventu- 
ally leading to a network splitting. Therefore, our model 
considers 2-dimcnsional lattices with different groups of 
nodes which randomly alternate between two possible 
states, i.e. they act as sources or sinks respectively. 
These state transitions induce time-varying stochastic 
flows on every link. Once reaching its capacity, a link 
fails with a time delay due to inertia effects. 

The motivation for this dynamic flow model was the 
large-scale integration of renewable intermittent energy 
sources (e.g. wind power, photovoltaic systems) into the 
electric power grid. This implies a higher ratio of non- 
dispatchable generation which, in turn, leads to less pre- 
dictable and more fluctuating flows on the network. Con- 
sequently, the anticipation of undesired situations such as 
cascading transmission line overloads leading to a net- 



work breakdown becomes highly complicated (Toj . In 
such a future infrastructure layout the network merely 
serves as a backbone for the redistribution of power 
from regions of energy surplus to regions with net power 
consumption. As detailed modeling and simulation ap- 
proaches become limited due to the increased complexity 
of electric power systems with large share of renewables, 
we opted for a minimalistic approach in order to under- 
stand the fundamental physics governing the dynamic 
behavior leading to a network splitting. Past experience 
has shown that such a splitting potentially results in a 
wide-area blackout with severe social and economic con- 
sequences fTi"| . 

Questions to be tackled are: What is the relation be- 
tween the stochastic behavior of the nodes and the emerg- 
ing flow patterns on the network? How are these flow 
patterns affected by different groupings of the nodes? 
What, in turn, is the impact of these flow patterns on 
the probability of a network splitting? How do inertia 
effects influence the potential splitting process? 

Although the definition of our model is based on the 
specific properties of future energy networks, it is ex- 
pected to reflect basic features of other real-world net- 
worked systems, whose robustness is subject to stochastic 
intermittent flows. 



II. DYNAMIC MODEL AND SIMULATION 
PROCEDURE 

Our study system incorporates a model for the nodal 
state alternation, a flow model, a lattice layout model 
and a model for the cascading link outages. 



A. Stochastic nodal state alternation 

The two possible states between which all nodes can al- 
ternate assume current injections of P± = 1 (node state 
"up") when the node acts as a source or Pr = — 1 (node 
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state "down" ) when the node acts as a sink. This stochas- 
tic up-down-up cycle assumes for every node i constant 
transition rates Xi and fii respectively. Hence, this alter- 
nating process is characterized by the cumulative distri- 
bution functions of the up-state and down-state times, 
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where t u and td are the time spans measured from the 
moment of entering the up-state and down-state respec- 
tively. The state transition frequency of every node is 
calculated by 



fi 
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(2) 



and corresponds to the average number of up-down-up 
cycles per time unit. For simplicity we assign to every 
node i the same transition rates Xi = X and fii = [/,, im- 
plying the same transition frequency fi = /. Moreover, 
the ratio is kept constant at X/ fi = 1 in order to assign 
the same probabilities to both possible states. 



B. Flow model 

We model the flows on the network by applying an elec- 
trical direct current model based on Ohm's law. Thereby, 
the linear relation between the nodal current injections 
Pi and the voltages Vi can be put into matrix form 



P = BV. 



(3) 

« - ■ a and 

is the resistance of each link 



The conductance matrix B has elements E>a = — 

R n = Ejen, r T^ where 

and fi, is the set of all the directly connected nodes 
to i. By assuming for simplicity that r,j = 1 for all links, 
the flow on a link is given by 



Pii=Vi 



Vi. 



(4) 



The sum of all the current injections at a given time in- 
stant is not necessarily equal to zero due to the stochastic 
nature of the up-down-up cycle. In order to satisfy the 
balance condition ^ Pj = at all times, a lack or sur- 
plus of the total current injections within the network is 
compensated by an additional, equally distributed injec- 
tion ±\(J2i Pi)/N\ at every node. Nevertheless, the sat- 
isfaction of the balance condition implies that the rows 
of B are linearly dependent. To make Eq. Q uniquely 
solvable, one of the equations in the system is removed 
and the node associated with that row is chosen as the 
voltage reference V re f = 0. 



C. Lattice layout and node grouping 



and L = 2N links with periodic (or "wrap-around") 
boundary conditions. In this way every node is directly 
connected to 4 neighbors, thus different conditions for 
boundary nodes are avoided. Furthermore, we partition 
the lattice into several square groups, each containing an 
equal number of nodes [Fig. [1] (a)]. All the nodes in a 
given group are in the same state at all times and al- 
ternate states simultaneously. We denote the grouping 
factor G as the number of groups in the network, thus 
G = N represents total stochastic independence between 
all nodes. As depicted in Fig. [1] (b)-(e), an increased / 
results in a higher fluctuation frequency of the flows. By 
further varying the grouping factor G, a broad spectrum 
of different stochastic flow patterns can be reproduced. 
A high value of G is leading to more smooth flow time- 
series, while a small value implies a strong fluctuation 
around the mean value. 



D. Link outage model 

In order to incorporate inertia effects in our model, the 
link outage mechanism is based on the concept that the 
flow Pij(t) determines the "temperature" on the 

link according to 



dT. l3 {t) 
dt 
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The link fails if Ty (t) reaches its capacity Tfj. In order 
to simplify Eq. ((SJ we set = 1. The parameter 
represents the characteristic time (inertia) constant. 

As an example, such an inertia is present in electric 
power grids where the power flows might heat the trans- 
mission lines up to a maximum allowable temperature. 



E. Simulation procedure 

With respect to the implementation we opted for a 
discrete-event based approach. This allows describing the 
time evolution of the nodal states and the resulting flows, 
as well as of the link outages and the resulting lattice 
status. By means of extensive Monte Carlo simulations 
we estimated the expected time until the splitting of the 
lattice. The simulation procedure comprises the following 
steps: 

1. Construct the N x N lattice adjacency matrix 
A and the N x N conductance matrix B. For 
all the nodes i in a single group determine their 
equal output states Pj at t = by a single 
Bernoulli trial with probability p = 0.5. Set the 



simulation step to n . — 0, set t 



(0) 



and ini- 



tialize the temperature of each link to Tij(trfy) = 0. 



We embedded our model for the nodal behavior and 
the resulting flows in a regular square lattice of N nodes 



2. Calculate the flow P^ on each link by Eq. (0| 
after solving Eq. ([3]) for V. For all links 
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FIG. 1. (Color online) (a) Schematic plot of a 12 x 12 lattice (JV = 144) with grouping factor G = 9 (left) and G = 36 (right). 
Note that the coloring of the groups has been used only due to illustrative reasons, all groups are stochastically independent 
with each other. The total flow \Pij\ versus time in the left lattice is depicted in (b) for / = 0.01 and in (d) for / = 0.5. 

Similar for the right lattice in (c) and (e) respectively. Notice that an increased / results in a higher fluctuation frequency of 
the flows. By further varying the grouping factor G, a broad spectrum of different stochastic flow patterns can be reproduced. 
A high value of G is leading to more smooth flow time-series, while a small value implies a strong fluctuation around the mean 
value. 



determine the subsequent time step Ai'^^^ after 
which they fail. If Pij(t n ) > X^, this time span is 
given by 



-Tij In 



(0) 



For every link calculate the point in time when 



temp 



t 



(») 



it fails due to reaching as t v] 
At ter ?P-, and build the vector t temp with ele- 
ments t t ^ mp . Determine the time of the first 
link outage as t out > temp = mm[t temp }. Deter- 
mine for every node i of the network the point 
in time if when it changes its state. Then, 
the time of the first state change is given by 

t change,s = mm [^] with ^ = % ... ^] _ De _ 

tcrmine the time of the next simulation event as 

t next = min { t out,te m p ^change,s^ Increment t h e sim- 
ulation step to n = 1. 

3. Proceed the simulation to ti n ) = t next . Remove the 
failed link (if any) from the lattice and update A 



and B. Recalculate the output Pj of each node i 
based on Eq. Recalculate the flow Pij and the 
temperature Tij on each link The flow P^ 

remains constant at least until the next event. The 
temperature Tij(t(n)) is given by 



T ij(t(n)) - P ij(t(n-1)) 



-At, 



+T l3 (t(n-i))e ^ ( "\ 
where At {n) = t (n) - t (n _ iy 
4. For each node and link recalculate if and t 



temp 



(7) 



and 



update t change ' 8 and t oui - tem P respectively as de- 
scribed in Step 2. Determine the time of the next 
event t next . 

5. Check the connectivity of the lattice. If it remains 
connected, increment the simulation step n and go 
back to Step 3. Otherwise, stop the simulation. 
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FIG. 2. (Color online) Average flows (Pij) versus the grouping 
factor G, which are well fitted by a power law with character- 
istic exponents b, slightly increasing with the lattice size TV. 
The dotted lines serve as a guide to the eye. 
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a higher current exchange among the groups. Moreover, 
(Pij) decreases with the number of independently alter- 
nating groups. For a given lattice size TV, increasing G 
implies less nodes in the groups and thus less exchange 
among them. Seen from a different angle, a high value 
of G means that less nodes behave simultaneously in the 
same way, leading to a more local current exchange and 
less flows in the lattice. In contrast, a low G induces 
higher flows over longer distances. Interestingly, for a 
given TV the decrease of the average flow with G follows 
a power law 



(P^) ocG- b ,G>9. 



(9) 



The exponent b is rather small and slightly increasing 
with the size of the lattice. 

As shown in Fig. [3] the data in Fig. [5] collapse onto 
a single curve, if the average flows are scaled with y/~N. 
This result can be explained by the flow distribution on 
the lattice. The average flow is largely determined by the 
maximum flows which are encountered at the boundaries 
of the groups. Suppose two lattices with sizes TVi and 
N2 and same grouping factor G. Then the maximum 
possible flows induced by a (square) group on one of its 
boundary links, p™ ax and P™ ax , are approximately pro- 
portional to y/Ni/G and \J N2/G with the same factor 
respectively. Thus p^xjpmax _ X /N 1 /N 2 . 



G 



FIG. 3. (Color online) Collapse of all the average flow data 



shown in Fig. [2] by scaling as (Pij) /vTV, versus the grouping 
factor G. The collapsed data follow a power law with charac- 
teristic exponent —0.27. The dotted lines serve as a guide to 
the eye. 



III. NUMERICAL RESULTS 



A. Average flows 

To clarify the impact of different grouping factors G 
and lattice sizes TV on the resulting flow patterns, we es- 
timate the average flow per link (Pij) without considering 
the link outage model, 

(Pij) = j Km (\ fY,\Pii{t')\dt' 

^EE(i p ^(«-i))i A V))> ( g ) 



t tot L 
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where t tot = J2 n ^« ^ s the sampled overall time span. 
The average flow is independent of /, as by increasing 
the frequency the relative duration among all different 
nodal state combinations remains unchanged. Figure [5] 
shows the values of (Pij) versus the grouping factor G 
in lattices of different sizes N. The average flows in- 
crease with the lattice size because for a given G the 
number of nodes in a group increases with TV leading to 



B. Lattice splitting 

The robustness of a lattice is quantified by the ex- 
pected time (t spin) when a splitting occurs and the lattice 
breaks into two parts [l2j . This time span can be inter- 
preted as the life expectancy of the lattice and depends 
on the capacity T?- of each link (i,j). To simplify mat- 
ters, we assign the same value T?- = T c and = r to all 
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FIG. 4. (Color online) Relationship between the link capacity 
T c and the expected time until the splitting, (t ap iu), of a 
24 x 24 lattice (TV = 576) for two different grouping factors G 
and state transition frequencies /. The inertia constant is set 
to r = 1. The error bars indicate the 95% confidence interval. 
The dotted lines serve as a guide to the eye. 
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FIG. 5. (Color online) Expected time {t sp ut) until the lattice splitting versus the state transition frequency /. If not stated 
otherwise, the lattices have size N = 576. The error bars indicate the 95% confidence interval, (a) Effect of different grouping 
factors G. The parameters of the dynamic model are set to r = 1 and T c — 4. (b) Collapse of the splitting time data in lattices 
of different sizes N with G = 36, by adjusting the link capacities according to T c — a^/N with a — 1/6. The inertia constant 
is set to t = 1. The average splitting times without adjusting T c are depicted in the inset for N = 576 and iV = 900. The 
corresponding values for N — 144 are omitted as the high splitting times induce prohibiting simulation run-times, (c) Collapse 
of (tsput) with t = 1 due to adjusting T c according to the grouping factor G. The inset shows the chosen value of T c for each 
value of G, being fitted by a power law. (d) Effect of the inertia constant r on the splitting times with G = 36 and T c — 4. 



links. Figure 2] shows the behavior of (t sp nt) versus an in- 
creasing value of T c for two different grouping factors G 
and state transition frequencies /. The expected time un- 
til the lattice splits increases exponentially with the link 
capacity. Hence, (t sp u t ) is highly sensitive with respect 
to small changes of T c . For a given value of T c , a larger 
grouping factor G leads to a significantly higher value of 
(t spin), as less flows are induced [Fig. [2]. The effect of 
varying the state transition frequency / is similarly large 
and is examined in more detail in Fig. [5] Starting with 
a low value, an increase of / leads to a shorter time span 
until the combined nodal states induce those minimum 
flows which are needed for the temperatures to reach 
the capacities T c [Fig. Q]. Consequently, as depicted 
in Fig. [5] (a)-(d), the splitting times (t sp i it ) are high at 
low values of / and become significantly decreased as the 



value of / is increasing. However, as / is exceeding a 
certain value, the splitting times start to increase again, 
and the lattice becomes more robust. This result can 
be explained by the inertia effects according to Eq. ([5]). 
While the flows on the lattice reach more often higher 
absolute values [Fig. [Tj, the average residence times of 
the underlying nodal states begin to fall below the mini- 
mum time needed to heat the links up to their capacity 

With a small number of groups the combined output is 
more fluctuating between the extreme values [Fig. (TJ left, 
compared to Fig. [1] right)]. This, in turn, is leading to a 
higher probability to encounter high flows and high link 
temperatures in a given time span. The splitting times 
(t S piit) thus are significantly shorter, as depicted in Fig. 
El (a). 
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By considering the scaling behavior of the average flows 
[Fig. [3J, the values of (t sp u t ) collapse for lattices with 
different sizes, but equal model parameters otherwise, if 
the link capacities are set as T c oc y/N. This result is 
demonstrated in Fig. [5](b) for three different lattice sizes. 
Notice that for a given value of / the average splitting 
times without adjusting T c differ by several decades [Fig. 
[5j inset]. 

As shown in Fig. [5] (c) the link capacities T c can be 
adjusted in such a way, that the splitting times in lat- 
tices with equal N but different grouping factors G over- 
lap for a wide range of the state transition frequency 
/. The adjusted capacities can be fitted by a power law 
with characteristic exponent —0.23 [Fig. 0(c), inset], be- 
ing remarkably close to the characteristic exponent b of 
Eq. ©. 

The effect of the inertia is shown in Fig. [5] (d) by 
varying the inertia constant r. Without any inertia, i.e. 
r = implying Zy(i) = P i:j (t) [Eq. ©], the splitting 
time declines with slope 1//. In average, the number of 
state transition events increases linearly with / in a given 
time span. This, in turn, decreases the average time until 
a maximum allowable flow Pjj = T c on a link is 
reached, in an inversely proportional manner. However, 
for t > a minimum average splitting time arises for a 
roughly estimated value of / « 0.05/t. 

In order to quantify the damage after the splitting, we 
follow [l3[ and measure the average relative size of the 
larger of the two remaining connected components 



(C) = (N') /N, 



(10) 



where (N') denotes the average number of nodes in the 
larger connected component. Figure [6] shows the average 
relative size of the larger connected component (C) for 
different grouping factors G and link capacities T c ver- 
sus the state change frequency /. The size of the larger 
connected component is rather independent of the flow 
fluctuation frequency as determined by /. While keeping 
the same link capacities [Fig. [51 values for T c = 4] there 
is no clear indication with regard to the dependence of 
(C) on the grouping factor G. For a given G, decreas- 
ing the link capacities T c slightly increases the value of 
(C) . For a smaller value of T c lower flows are sufficient to 
overload the links and cascades may develop in smaller 
regions. Hence, the failing links envelop a lower num- 
ber of nodes eventually breaking away from the lattice, 
implying a larger size of the remaining connected compo- 
nent after the splitting. However, for the chosen values 
of T c the average relative size remains approximately in 
the range 0.65 < (C) < 0.75. 



IV. SUMMARY AND CONCLUSIONS 

To sum up, in this paper we have introduced a min- 
imal model for stochastic intermittent flows on lattices. 
These flows might induce cascading link overloads even- 
tually leading to a lattice splitting into two parts. In or- 
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FIG. 6. (Color online) Average size of the larger connected 
component (C) after the splitting into two parts versus the 
state transition frequency / for various grouping factors G 
and link capacities T c . The lattices have size N = 576. The 
error bars indicate the 95% confidence interval. The dotted 
lines serve as a guide for the eye. 



der to better represent real systems we implied an inertia 
in such a way that a link does not fail immediately but 
rather delayed when it becomes overloaded. By extensive 
Monte Carlo simulations we revealed how the time until 
such a splitting occurs depends on different flow patterns. 
With an increasing number of (stochastically) indepen- 
dent nodes the average flows decrease slowly, following a 
power law. With regard to the robustness of the lattices, 
a high sensitivity of the splitting time to the link capaci- 
ties is observed. Increasing the flow fluctuation frequency 
(as determined by the nodal state alternation) decreases 
this time span until reaching a minimum, after which it 
rises again meaning a higher "life expectancy" of the lat- 
tice. Generally, both a higher stochastic independence 
among the nodes (i.e. more groups of simultaneously al- 
ternating nodes) and a smaller size of the lattice imply 
higher splitting times. However, these time spans seem 
to coincide by adjusting the link capacities according to 
a power law with respect to the node grouping, and ac- 
cording to the square-root of the lattice size respectively. 
Furthermore, we have shown that the effect of the inertia 
is significant. Its absence implies a monotonic decrease 
of the splitting times, while introducing it results in re- 
markably higher values for higher inertia constants. As 
an indication of the damage after the splitting, the rela- 
tive size of the larger connected component seems to be 
independent of the flow fluctuation frequency but sligthly 
decreases with the link capacity. 

We conclude with some thoughts on the implications 
of these results for future energy networks, being charac- 
terized by a large share of renewable intermittent power 
sources. The more distributed the power sources are (be- 
ing equivalent to more groups in our model), the lower 
the flows exchanged over the power grid [Fig. [2] and Fig. 
[5] (a)]. However, even in a highly distributed system, a 
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considerable transmission capacity is still needed to keep 
the system at the desired level of security [Fig. [5] (c)]. 
Increasing the size of the grid can be expected leading 
to a disproportionally small increase of the flows [Fig. 
|3] ■ Restricting the capacities of the transmission lines 
or, equivalently, operating the system closer to its secu- 
rity margins might reduce the robustness of the network 
against cascading failures drastically [Fig. [3]. The inertia 
as induced by the heating of the transmission lines which 
might fail when reaching a maximum allowable tempera- 
ture, potentially increases the robustness of power grids 
with large share of renewables [Fig. [5] (d)]. The same 
effect can be even exploited for increasing existing trans- 
mission line capacities, thus improving the economic per- 
formance of the system (3] ■ If the grid breaks apart as 
a result of cascading line failures, the sizes of the two 
formed islands can be expected to be largely indepen- 
dent of the flow fluctuation frequencies. Nevertheless, 
they seem to slightly become more asymmetric with de- 



creasing power transfer capacities thus leaving a larger 
remaining connected component of the network [Fig. [B]. 

Our model provides insights into the underlying 
physics of networks subject to stochastic flows. There- 
fore, we believe that besides future energy networks, po- 
tential applications could be investigated on other real- 
world systems such as traffic networks. 
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